setwd("06.phenotype/")
options(stringsAsFactors = F)
library(dplyr)
library(ggplot2)
library(reshape2)
raw <- read.delim("quality",header = T,
                  stringsAsFactors = F,
                  sep = "\t")
raw$MIC = as.numeric(gsub(pattern = "-",replacement = 0, raw$MIC))
str(raw)
## 按组合分开
BtJ <- filter(raw, Sample == "Bt" | Sample == "BtJ" | Sample == "J")
BtSY <- filter(raw, Sample == "Bt" | Sample == "BtSY" | Sample == "SY")
Z12B1 <- filter(raw, Sample == "Z12" | Sample == "Z12B1" | Sample == "B1")
Z41S <- filter(raw, Sample == "Z41" | Sample == "Z41S" | Sample == "S")
B985S <- filter(raw, Sample == "B985" | Sample == "B985S" | Sample == "S")
L28ZA <- filter(raw, Sample == "L28ZA" | Sample == "L28" | Sample == "ZA")
Z12J <- filter(raw, Sample == "Z12" | Sample == "Z12J" | Sample == "J")
## 放到一个list中
filename = c("BtJ","BtSY","Z12B1","Z41S","B985S","L28ZA","Z12J")
raw.data <- list(BtJ,BtSY,Z12B1,Z41S,B985S,L28ZA,Z12J)
clean.data <- list()
head(BtJ)
for (i in 1:length(filename)) {
  raw.data[[i]]$Type = if_else(raw.data[[i]]$Sample == "Bt" | raw.data[[i]]$Sample == "Z12" | raw.data[[i]]$Sample == "Z41" | raw.data[[i]]$Sample == "L28" | raw.data[[i]]$Sample == "B985","Female",
                               if_else(raw.data[[i]]$Sample == "BtJ" | raw.data[[i]]$Sample == "BtSY" | raw.data[[i]]$Sample == "Z12J" | raw.data[[i]]$Sample == "Z12B1" | raw.data[[i]]$Sample == "Z41S" | raw.data[[i]]$Sample == "B985S" | raw.data[[i]]$Sample == "Z41S" | raw.data[[i]]$Sample == "L28ZA","Hybrid",
                                       if_else(raw.data[[i]]$Sample == "J" | raw.data[[i]]$Sample == "SY" | raw.data[[i]]$Sample == "ZA" | raw.data[[i]]$Sample == "B1" | raw.data[[i]]$Sample == "S","Male","0")))
  raw.data[[i]]$Sample = rep(filename[i],nrow(raw.data[[i]]))
  
}
## 合并
data.merge <- bind_rows(raw.data)
head(data.merge)
write.table(data.merge, "~/06.phenotype/quality_Datamelt.xls", sep = "\t", row.names = F, quote = F)
Trait = colnames(raw.data[[1]])[2:6]
label = c("cm","%","","cN/Tex","%")
for (i in 1:length(Trait)) {
  p = ggplot(data.merge, aes(x = Sample, y = data.merge[,i+1], fill = Type)) + 
    geom_boxplot(position = position_dodge(0.8)) +
    geom_dotplot(binaxis = 'y', stackdir = "center",,dotsize = 0.2,position = position_dodge(0.8)) +
    labs(title = Trait[i], x = "Cross combination", y = label[i]) +
    theme_bw()
  ggsave(plot = p, filename = paste(Trait[i],"boxplot.pdf",sep = "-"),width = 15, height = 10)
}

head(BtJ)

Heterosis <- function(data, P1, P2, F1){
  Pmean <- apply(filter(.data = data, Sample == P1 | Sample == P2)[,-1], 2, mean)
  Fmean <- apply(filter(.data = data, Sample == F1)[,-1], 2, mean)
  Hbp = 100*(Fmean - Pmean)/Pmean
  P1.value = apply(filter(data, Sample == P1)[,-1],2,mean)
  P2.value = apply(filter(data, Sample == P2)[,-1],2,mean)
  F1.value = apply(filter(data, Sample == F1)[,-1],2,mean)
  Pplus = if_else(P1.value > P2.value, P1.value,P2.value)
  Sbp = 100*(F1.value-Pplus)/Pplus
  x <- as.data.frame(rbind(Hbp,Sbp))
  x$Type = c("MPV","OPV")
  z = melt(data = x,variable.names="Type",value.name = "Num")
  write.table(z,file = paste(F1,"MPHandOPH.xls",sep = ""),
              row.names = F,sep = "\t",quote = F)
  p = ggplot(data = z, aes(x = z[,2], y = Num, fill = Type)) +
    geom_bar(stat ="identity",width = 0.6,position = "dodge") +
    labs(title = "Mid-parent heterosis and Over-parent heterosis", x = F1, y = "%") +
    theme_bw()
  ggsave(p,filename = paste(F1,"MPHandOPH.pdf",sep = "-"),width = 15,height = 10)
  assign("z",value = z, envir = globalenv())
}
filename = c("BtJ","BtSY","Z12B1","Z41S","B985S","L28ZA","Z12J")
raw.data <- list(BtJ,BtSY,Z12B1,Z41S,B985S,L28ZA,Z12J)
P1 = c("Bt","Bt","Z12","Z41","B985","L28","Z12")
P2 = c("J","SY","B1","S","S","ZA","J")
F1 = c("BtJ","BtSY","Z12B1","Z41S","B985S","L28ZA","Z12J")

het <- list()

for (i in 1:length(F1)) {
  Heterosis(data = raw.data[[i]],
            P1 = P1[i],
            P2 = P2[i],
            F1 = F1[i])
  het[[i]] <- z
  het[[i]]$Var <- F1[i]
}

het1 <- bind_rows(het)
het1
MPV <- filter(het1, Type == "MPV")[,-1]
OPV <- filter(het1, Type == "OPV")[,-1]
head(MPV)

write.table(MPV,file = "MPV30Ind.xls",row.names = F, sep = "\t", quote = F)
write.table(OPV,file = "OPV30Ind.xls", row.names = F, sep = "\t", quote = F)


